
## setwd("C:/Users/Jon/Documents/aid/aidnicolas")

library(riv)
library(foreign)

set.seed(101, kind = "default", normal.kind = NULL)

## work out threshold (for static model)

## dthresh <- sqrt(qchisq(.50, df=990))
## dthresh <- sqrt(qchisq(.99, df=3))
dthresh <- 20

## NB note rename csv file so no spaces in filename

nicdata3 <- read.table("StaticFSandSS.csv",sep=",",header=T)
attach(nicdata3)

Y2 <- as.matrix(nicdata3[,7]) 
Xend2 <- as.matrix(nicdata3[,8]) 
Zinst2 <- as.matrix(nicdata3[,9])
Xex2 <- rep(1,length(Y2))

civcce <- riv(Y2, Xend2, NULL, Zinst2, method="classical")
crob <- riv(Y2, Xend2, NULL, Zinst2, method='S-est')
dcrob <- sqrt(crob$MD)
cdrop <- country[dcrob>dthresh]
wmat <- cbind(obsnum,crob$weight)
distmat <- cbind(obsnum,ccode,as.vector(dcrob))

## Government consumption

Y <- as.matrix(nicdata3[,10]) 
Xend <- as.matrix(nicdata3[,11]) 
Zinst <- as.matrix(nicdata3[,12])

givcce <- riv(Y, Xend, NULL, Zinst, method="classical")
grob <- riv(Y, Xend, NULL, Zinst, method="S-est")
dgrob <- sqrt(grob$MD)
gdrop <- country[dgrob>dthresh]
wmat <- cbind(wmat,grob$weight)
distmat <- cbind(distmat,as.vector(dgrob))

## Total consumption

Y <- as.matrix(nicdata3[,13]) 
Xend <- as.matrix(nicdata3[,14]) 
Zinst <- as.matrix(nicdata3[,15])

tivcce <- riv(Y, Xend, NULL, Zinst, method="classical")
trob <- riv(Y, Xend, NULL, Zinst, method="S-est")
dtrob <- sqrt(trob$MD)
tdrop <- country[dtrob>dthresh]
wmat <- cbind(wmat,trob$weight)
distmat <- cbind(distmat,as.vector(dtrob))

## Gross capital formation

Y <- as.matrix(nicdata3[,16]) 
Xend <- as.matrix(nicdata3[,17]) 
Zinst <- as.matrix(nicdata3[,18])

iivcce <- riv(Y, Xend, Xex=NULL, Zinst, method="classical")
irob <- riv(Y, Xend, Xex=NULL, Zinst, method="S-est")
dirob <- sqrt(irob$MD)
idrop <- country[dirob>dthresh]
wmat <- cbind(wmat,irob$weight)
distmat <- cbind(distmat,as.vector(dirob))

## Exports

Y <- as.matrix(nicdata3[,19]) 
Xend <- as.matrix(nicdata3[,20]) 
Zinst <- as.matrix(nicdata3[,21])

xivcce <- riv(Y, Xend, Xex=NULL, Zinst, method="classical")
xrob <- riv(Y, Xend, Xex=NULL, Zinst, method="S-est")
dxrob <- sqrt(xrob$MD)
xdrop <- country[dxrob>dthresh]
wmat <- cbind(wmat,xrob$weight)
distmat <- cbind(distmat,as.vector(dxrob))

## Imports

Y <- as.matrix(nicdata3[,22]) 
Xend <- as.matrix(nicdata3[,23]) 
Zinst <- as.matrix(nicdata3[,24])

mivcce <- riv(Y, Xend, Xex=NULL, Zinst, method="classical")
mrob <- riv(Y, Xend, Xex=NULL, Zinst, method="S-est")
dmrob <- sqrt(mrob$MD)
mdrop <- country[dmrob>dthresh]
wmat <- cbind(wmat,mrob$weight)
distmat <- cbind(distmat,as.vector(dmrob))

## Net imports

Y <- as.matrix(nicdata3[,25]) 
Xend <- as.matrix(nicdata3[,26]) 
Zinst <- as.matrix(nicdata3[,27])

mxivcce <- riv(Y, Xend, Xex=NULL, Zinst, method="classical")
mxrob <- riv(Y, Xend, Xex=NULL, Zinst, method="S-est")
dmxrob <- sqrt(mxrob$MD)
mxdrop <- country[dmxrob>dthresh]
wmat <- cbind(wmat,mxrob$weight)
distmat <- cbind(distmat,as.vector(dmxrob))

dropallidx <- ((dcrob>dthresh)|(dgrob>dthresh)|(dtrob>dthresh)|(dirob>dthresh)|(dxrob>dthresh)|(dmrob>dthresh)|(dmxrob>dthresh))
dropall <- country[dropallidx==1]
dropnum <- obsnum[dropallidx==1]
distord <- order(dcrob,decreasing = TRUE)
distmat <- as.data.frame(distmat[distord,])
distmat <- cbind(country[distord],distmat)

colnames(wmat) <- c("obsnum","cw","gw","tw","iw","xw","mw","mxw")
write.dta(as.data.frame(wmat), "weightsdata.dta")

## now save Mahalanobis distances, ordered by distances for static household consumption regression

dimnames(distmat)[[2]] <- c("country","obsnum","ccode","cd","gd","td","id","xd","md","mxd")
write.dta(as.data.frame(distmat), "distancesdata.dta")
## write.xls(distmat, "distancesdata.xls")
write.table(distmat, file = "distancesdata.csv", sep = ",", col.names = NA, qmethod = "double")

sink("rivoutput.txt")

cat("Consumption")
civcce$Summary.Table
crob$Summary.Table

cat("Govt consumption")
givcce$Summary.Table
grob$Summary.Table

cat("Total consumption")
tivcce$Summary.Table
trob$Summary.Table

cat("Investment")
iivcce$Summary.Table
irob$Summary.Table

cat("Exports")
xivcce$Summary.Table
xrob$Summary.Table

cat("Imports")
mivcce$Summary.Table
mrob$Summary.Table

cat("Net imports")
mxivcce$Summary.Table
mxrob$Summary.Table

dropall
dropnum

sink()

## experiments

Y <- as.matrix(nicdata3[,7]) 
const <- rep(1,nrow(Y))
X <- as.matrix(cbind(const,as.matrix(nicdata3[,8]))) 
Z <- as.matrix(cbind(const,as.matrix(nicdata3[,9])))

iv2sls <- solve(crossprod(Z,X),crossprod(Z,Y))

## now try weighted 2sls

sigma <- diag(crob$weight)
X2 <- sigma%*%X
Y2 <- sigma%*%Y
riv <- solve(crossprod(Z,X2),crossprod(Z,Y2))

Z2 <- crossprod(Z,sigma)
riv2 <- solve(Z2%*%X,Z2%*%Y)

## now try Gabriela Cohen-Freue's code

W <- diag(crob$weight)
GX <- as.matrix(nicdata3[,8])
GY <- as.vector(unlist(nicdata3[,7]))

X.til <- cbind(rep(1,length(GX)),GX)
Z.til <- cbind(rep(1,length(GX)),as.matrix(nicdata3[,9]))

zx.weight <- t(Z.til)%*%W%*%X.til
zx.weight <- solve(zx.weight)

zy.weight <- t(Z.til)%*%W%*%GY
rivwls <- zx.weight%*%zy.weight

detach("nicdata3")

## now try dynamic panel

nicdata4 <- read.table("DynamicFSandSS.csv",sep=",",header=T)
attach(nicdata4)

Y2 <- as.matrix(nicdata4[,7]) 
Xend2 <- as.matrix(nicdata4[,8]) 
Zinst2 <- as.matrix(nicdata4[,9])
LY2 <- as.matrix(nicdata4[,10])
Xex2 <- rep(1,length(Y2))

civcce <- riv(Y2, Xend2, LY2, Zinst2, method="classical")
crob <- riv(Y2, Xend2, LY2, Zinst2, method='S-est')
dcrob <- sqrt(crob$MD)
cdrop <- country[dcrob>dthresh]
wmat <- cbind(obsnum,crob$weight)
distmat <- cbind(obsnum,ccode,as.vector(dcrob))

## Government consumption

Y <- as.matrix(nicdata4[,11]) 
Xend <- as.matrix(nicdata4[,12]) 
Zinst <- as.matrix(nicdata4[,13])
LY2 <- as.matrix(nicdata4[,14])

givcce <- riv(Y, Xend, LY2, Zinst, method="classical")
grob <- riv(Y, Xend, LY2, Zinst, method="S-est")
dgrob <- sqrt(grob$MD)
gdrop <- country[dgrob>dthresh]
wmat <- cbind(wmat,grob$weight)
distmat <- cbind(distmat,as.vector(dgrob))

## Total consumption

Y <- as.matrix(nicdata4[,15]) 
Xend <- as.matrix(nicdata4[,16]) 
Zinst <- as.matrix(nicdata4[,17])
LY2 <- as.matrix(nicdata4[,18])

tivcce <- riv(Y, Xend, LY2, Zinst, method="classical")
trob <- riv(Y, Xend, LY2, Zinst, method="S-est")
dtrob <- sqrt(trob$MD)
tdrop <- country[dtrob>dthresh]
wmat <- cbind(wmat,trob$weight)
distmat <- cbind(distmat,as.vector(dtrob))

## Gross capital formation

Y <- as.matrix(nicdata4[,19]) 
Xend <- as.matrix(nicdata4[,20]) 
Zinst <- as.matrix(nicdata4[,21])
LY2 <- as.matrix(nicdata4[,22])

iivcce <- riv(Y, Xend, LY2, Zinst, method="classical")
irob <- riv(Y, Xend, LY2, Zinst, method="S-est")
dirob <- sqrt(irob$MD)
idrop <- country[dirob>dthresh]
wmat <- cbind(wmat,irob$weight)
distmat <- cbind(distmat,as.vector(dirob))

## Exports

Y <- as.matrix(nicdata4[,23]) 
Xend <- as.matrix(nicdata4[,24]) 
Zinst <- as.matrix(nicdata4[,25])
LY2 <- as.matrix(nicdata4[,26])

xivcce <- riv(Y, Xend, LY2, Zinst, method="classical")
xrob <- riv(Y, Xend, LY2, Zinst, method="S-est")
dxrob <- sqrt(xrob$MD)
xdrop <- country[dxrob>dthresh]
wmat <- cbind(wmat,xrob$weight)
distmat <- cbind(distmat,as.vector(dxrob))

## Imports

Y <- as.matrix(nicdata4[,27]) 
Xend <- as.matrix(nicdata4[,28]) 
Zinst <- as.matrix(nicdata4[,29])
LY2 <- as.matrix(nicdata4[,30])

mivcce <- riv(Y, Xend, LY2, Zinst, method="classical")
mrob <- riv(Y, Xend, LY2, Zinst, method="S-est")
dmrob <- sqrt(mrob$MD)
mdrop <- country[dmrob>dthresh]
wmat <- cbind(wmat,mrob$weight)
distmat <- cbind(distmat,as.vector(dmrob))

## Net imports

Y <- as.matrix(nicdata4[,31]) 
Xend <- as.matrix(nicdata4[,32]) 
Zinst <- as.matrix(nicdata4[,33])
LY2 <- as.matrix(nicdata4[,34])

mxivcce <- riv(Y, Xend, LY2, Zinst, method="classical")
mxrob <- riv(Y, Xend, LY2, Zinst, method="S-est")
dmxrob <- sqrt(mxrob$MD)
mxdrop <- country[dmxrob>dthresh]
wmat <- cbind(wmat,mxrob$weight)
distmat <- cbind(distmat,as.vector(dmxrob))

dropallidx <- ((dcrob>dthresh)|(dgrob>dthresh)|(dtrob>dthresh)|(dirob>dthresh)|(dxrob>dthresh)|(dmrob>dthresh)|(dmxrob>dthresh))
dropall <- country[dropallidx==1]
dropnum <- obsnum[dropallidx==1]
distord <- order(dcrob,decreasing = TRUE)
distmat <- as.data.frame(distmat[distord,])
distmat <- cbind(country[distord],distmat)

colnames(wmat) <- c("obsnum","cw","gw","tw","iw","xw","mw","mxw")
write.dta(as.data.frame(wmat), "dynamicweightsdata.dta")

## now save Mahalanobis distances, ordered by distances for dynamic household consumption regression

dimnames(distmat)[[2]] <- c("country","obsnum","ccode","cd","gd","td","id","xd","md","mxd")
write.dta(as.data.frame(distmat), "dynamicdistancesdata.dta")
## write.xls(distmat, "dynamicdistancesdata.xls")

write.table(distmat, file = "dynamicdistancesdata.csv", sep = ",", col.names = NA, qmethod = "double")

sink("rivoutputdynamic.txt")

cat("Consumption")
civcce$Summary.Table
crob$Summary.Table

cat("Govt consumption")
givcce$Summary.Table
grob$Summary.Table

cat("Total consumption")
tivcce$Summary.Table
trob$Summary.Table

cat("Investment")
iivcce$Summary.Table
irob$Summary.Table

cat("Exports")
xivcce$Summary.Table
xrob$Summary.Table

cat("Imports")
mivcce$Summary.Table
mrob$Summary.Table

cat("Net imports")
mxivcce$Summary.Table
mxrob$Summary.Table

dropall
dropnum

sink()
